Numerical modeling of the central spin problem using the spin coherent states 

P-representation 
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In this work, we consider decoherence of a central spin by a spin bath. In order to study the non- 
perturbative decoherence regimes, we develop an efficient mean-field-based method for modeling 
the spin-bath decoherence, based on the P-representation of the central spin density matrix. The 
method can be applied to longitudinal and transversal relaxation at different external fields. In 
particular, by modeling large-size quantum systems (up to 16000 bath spins), we make controlled 
predictions for the slow long-time decoherence of the central spin. 
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PACS numbers: 75.40.Mg, 73.21. La, 75.10.Jm 

Detailed understanding of the decoherence of various 
quantum systems is important for many areas, from 
quantum optics and solid-state physics to the quantum 
information processing, where decoherence constitutes a 
major obstacle for building a practical quantum com- 
puter. For example, in a quantum dot-based architec- 
ture, the quantum bit is represented by a spin of a single 
electron (central spin) placed in a quantum dot (QD). 
Due to interaction with the bath of nuclear spins in a 
QD, the electron spin quickly "looses memory" about its 
initial orientation and can not be used for computation. 
Experimental studies of this process has become possible 
very recently 0, Q , and detailed theoretical understand- 
ing of the experimental results is timely and important. 
Moreover, the problem of a central spin coupled to the 
spin bath 3] (the "quantum central spin problem") has 
recently arised in other contexts (decoherence in mag- 
netic molecules, dynamics of the cold fermions pairing), 
and has attracted much attention. 

Decoherence is a complex quantum many-body phe- 
nomenon, and satisfactory solutions can be obtained only 
in very special cases 0, 0- Perturbation theory can 
be successfully applied in the case of a strong magnetic 
field or polarized nuclear spin bath (which produces a 
strong Overhauser field acting on the central spin) 0, Q • 
But for the experimentally important non-perturbative 
regimes, no well-justified method, numerical or analyti- 
cal, has been suggested yet. The static approximation 
\li L3 lil , which simply neglects the dynamics of the bath 
spins, works only at short times, while the interesting 
long-time dynamics of the central spin remains an open 
problem, and the simulations on moderate-sized systems 
(~ 20 spins) do not give conclusive results. 

In this work, we present a novel approach to the quan- 
tum central spin problem based on the time-dependent 
mean field (TDMF) theory. It has been pointed out pre- 
viously [lfl Hll | that the mean field approach should be 
adequate, since the central spin interacts with a large 
number N of the bath spins (loosely speaking, the num- 



ber of the "nearest neighbors" for the central spin is 
large, i.e. the coupling of the central spin with the bath 
has an infinite range). However, the standard mean field 
approach |l2| gives a clearly wrong answer (see below). 
We use the spin coherent state P-representation 0] to 
modify the standard TDMF, and present an efficient ap- 
proach, which gives excellent agreement with the exact 
solution of the many-spin quantum problem already for 
rather small systems (up to 20 spins). By applying this 
method to large-scale problems, we study the interesting 
long-time dynamics of the central spin. Moreover, the 
P-representation approach allows us to understand why 
the properly corrected TDMF theory works for large N, 
and what the limitation of the method are. 

It is interesting to point out that the spin coher- 
ent states are traditionally used in the spin path inte- 
grals and in the semiclassical approximation for quan- 
tum spins. The powerful methods based on P- and Q- 
representations, so useful in quantum optics, have not 
been widely applied to the spin systems studies. We ex- 
pect that this work will help developing novel coherent- 
state approaches to studying spin systems. Another in- 
teresting point concerns the basic ideas of the mean field 
theory. The standard mean field theory approximates, 
in the optimal way, the exact many-spin wavefunction as 
a product of single-spin wavefunctions. However, when 
studying decoherence, we are only interested in the state 
of the central spin, not the whole many-body state. It 
is interesting that there is a justified modification of 
the standard TDMF (presented below) which provides a 
much better approximation for the relevant observables 
(the state of the central spin) at the expense of irrelevant 
information (the state of the bath). 

The electron spin in a QD interacting with the bath of 
N nuclear spins can be described by the Hamiltonian 
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where S = (Sq,Sq,Sq) are the operators of the elec- 
tron spin, Sfc are the operators of the bath spins, and 
A k = (8ir /3)5*^sffnA'n| x I'( x fc)| 2 is the contact hypcrfinc 
coupling which is determined by the electron density 
|<3>(xfc)| 2 at the site x& of the k-th nuclear spin and by 
the Lande factors of the electron g* and of the nuclei g n . 
The first term in Eq.QJdescribes the Zeeman energy of the 
electron spin in the external magnetic field Bq, and the 
second term represents the contact hyperfine coupling. 
The omitted terms, such as the Zeeman energy of the 
nuclear spins and the non-isotropic part of the hyperfine 
coupling, are very small and can be neglected for a wide 
range of experimental situations. Eq. ^ is the standard 
Hamiltonian of the quantum central spin problem jjj. 

We are interested in the dynamics of the central spin, 
i.e. in the dynamics of So(i) = Trp(i)So, where pit) is 
the density matrix of the total system (central spin plus 
the bath). Although the quantum central spin problem 
is integrable, the formal solution Q is very complex, and, 
to our knowledge, it has not been used for actual calcula- 
tions of so(t), neither analytically, nor numerically. Effi- 
cient approximate approaches are needed, and the TDMF 
theory is a natural first step. Within the mean-field ap- 
proach, we approximate the wavefunction of the total 
system as a product |\I>) = \i[)q) ®k=\ iV'fc) of the single- 
spin wavefunctions The TDMF equations of motion 
for \ipj) are obtained by substituting this ansatz into the 
Dirac's functional D = J dt(^\i(d^ /dt) - W$>), and re- 
quiring that SD = with respect to small variations of all 
\ipk}- The resulting equations of motion can be presented 
in a simple form as 

a, = [hjxsj] (j=0...N) (2) 
h = H e z + ^2A k s k , h k =A k s Q (k = 1 . . . N)(3) 

k 

where Sj = Trp(t)Sj, and [hj x Sj] is the vector prod- 
uct of hj and Sj. The physical meaning of these equa- 
tions is simple: every j-th spin precesses in its own time- 
dependent effective field hj given by Eq. [3] However, 
TDMF theory gives a very bad approximation to the ex- 
act solution of the problem. It can be seen, for example, 
from comparison of the standard TDMF theory with the 
exact numerical solutions for several distributions of 
A k (see Fig.EJ: the disagreement is significant. 

In order to construct a working approximation based 
on the TDMF, let us consider the P-representation of 
the system's density matrix in the basis of spin coherent 
states 0, . The spin coherent state for spin J is de- 
fined as = A^i J(2 J)!/(m!( J-m)!)] Vy™| J_ m ), 
where W = (1 + |/i| 2 ) J is the normalization constant. 
For a spin 1/2, the coherent state has a simple form 
\fi) = cos (6»/2)| t) + sin(6»/2)e^| |), where we used the 
parametrization p = tan(6*/2)e^. The basis of coher- 
ent states is overcomplete, and by chosing an appro- 
priate real- valued function A(0,<f>), any hermitian op- 



erator A can be represented in a diagonal form: A = 
f„A(9, <f))\fi)(p,\ $va,9d9d(j>, where the integration is per- 
formed over the sphere. Note that A(9, <f>) ^ (9, 4>\A\9, <j>). 
Moreover, the function A(9, <p) is not uniquely deter- 
mined; if we add to it any linear combination of the 
spherical harmonics Y™{0, 4>) of the order I > 2 J, then 
the value of the integral remains the same, and the new 
function would define the same operator A. Obviously, 
the many-spin density matrix p can also be written in 
the diagonal form: 




(4) 

with a real- valued function p({0j, 4>j\, i), where {9j,(j>j} 
is the set of all 6q, . . . On and <f>o, . . . (pN- This representa- 
tion for the density matrix is called the P-representation 
fliij l . Note that the operator part of the expression (@J has 
a mean-field form, i.e. it is a product of single-spin den- 
sity matrices \pj){pj\. In P-representation, the quantum- 
mechanical average x — Tr(p(t)X) of any observable X 
can be calculated by a simple formula 

. N N 

x = /p({^> J },t)(g)(^|X| Mj )I]sin^^d0 J , (5) 

3=0 3=0 

Our goal is to model the evolution of the function 
p({9j,4>j},t). However, the direct solution of the com- 
plex partial differential equation for p{{Qj,<f>j},t) is im- 
possible for a large (hundreds or more) number of spins, 
and we use a different approach. We note that if 
p({0j,4>j},t) > then this function can be interpreted 
as a probability for the system to be in the product 
state \ty) = ®jLo I^j)> an< ^ we nee d to simulate the dy- 
namics of the probability distribution p({9j,(j>j},t). To 
do that, we initially generate many realizations of the 
random vector (6q , . . . 6^' , </)$ , . . . <j)ffl ) distributed 
according to the probability distribution p({9j,4>j},0) 
(the index m — 1, . . . M enumerates the different 
realizations). Then we propagate every initial vec- 
tor (0Q™ 1 ', . . . 9f^\<^Q, . . . 4>^ ) in time, so that af- 
ter a lapse of time t, the initial vector evolves into 
a vector (O m) (i), • ■ ■ (t), $ m) (i), . . . d>^(t)). If 
the equations of motion for all the variables Q^\t) 

and ^ n \t) are chosen correctly, then the func- 
tion p({8j,(j)j},t) — p(Oj(t),$j(t)), and the value 
x = Tr(p{t)X) of any observable X can be cal- 
culated as an average over all realizations: x = 

To implement this approach, we need to determine 
the equations of motion for Qj(t), $j(t) which would 
produce a good approximation for p({9j,<j>j},t). As a 
first step, let us find the exact equations of motion for 
p({9j, 4>j}, t). For simplicity, let us study one term in the 
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central-spin Hamiltonian £[J, i.e. we consider two spins 
1/2 (the central spin and the A:-th bath spin) coupled 
by the isotropic Heisenberg interaction Hk = AkS Sk- 
The most general form for the two-spin density matrix is 
p = woololfe + w 0a l cr% + wpocrQ lfc + w\ v a^a^, where 
a = x,y,z (and similarly for other Greek indices), and 
erg and erf denote the Pauli matrices for the 0-th and the 
k-ih spin, respectively. Here and below, we assume sum- 
mation over the repeating indices. From the von Neu- 
mann's equation p(t) — i[p(t),H], we obtain woo(t) = 
(which expresses that Tr p(t) = 1), and 

w a p{t) = — e a pj[wjo(t) - w 0j (t)}, (6) 

where e Q( a 7 is a completely antisymmetric unity tensor 
(permutation symbol) . These equations of motion deter- 
mine the dynamics of p({9q, <pQ, Ok, 4>k\, t). From the P- 
representation |0J it follows that p({9o,cf>o,9k,(f>k},t) — 
Poo(t) +Po a (t)c^ +P/3o(*)cq + P\v(t)c^c u k , where p Q0 = 
(l/47T 2 )w o, Poa = (3/47r 2 )wo Q , p a0 = (3/47r 2 )w a0 , and 
Pa/3 — (9/47r 2 )w Qi (3. Here, we used the shorthand nota- 
tions Cq — sin 9q cos 4>o, Cq = sin do sin fa, Cq = cos 9q 
(and similarly for c|, ct, c|). The spherical harmon- 
ics of the order higher than one can also be added to 
p({0q, (f)Q, 9k, 4>k}), but they do not change the density 
matrix p, and therefore are not physically significant. 

The mean-field structure of the P-representation for 
the density matrix (Eq. |3J| suggests that the equations of 
motion for {Qo(t), $>o{t), ®k(t)} should also have 

a mean-field form corresponding to Eq. [3 but the local 
fields should be re-defined to provide optimal approx- 
imation for so(f). For simplicity, we omit the discus- 
sion of the general form for ho(t) and hk(t), and proceed 
to the answer. We introduce the shorthand notations 
Cq = sin6 cos $ , Cq = sin6osin$ 0l Cq = cos6 
(and similarly for C£ , C\, Cf ), and postulate the follow- 
ing equations of motion: 

C = g x [C fc x C ] , C k = -92 [C fc x C ] (7) 

(cf. Eq. EJ), where g\ and gi are to be determined. By 
substituting these equations into the probability distri- 
bution p({& (t), $o(t), e fc (i), **(*)}) - -Poo + P 0a C^ + 
P/3oCq + P\ v CqC^., and using the P-representation 
we obtain the following equations of motion for the den- 
sity matrix p: 

Wo 7 (t) = -g2ea/37«Vj(i), w y0 (t) = gieufjyWapit), 
w a p{t) = (l/3)e Q/ 3 7 [g 2 wo 7 (i) - giw l0 (t)}, (8) 

(cf. Eqs.EJ, and, trivially w 00 (t) = 0. The Eqs. © and0 
are incompatible for any g\ and g2 (i.e., TDMF is never 
exact). However, we are interested only in w-yo(t), since 




Time Time 

FIG. 1: Longitudinal Sg(t) (a) and transversal s%{t) (b) re- 
laxation of the central spin 1/2 coupled to a bath of N = 21 
spins 1/2. Couplings Ak are randomly distributed between 
and 1.0, external fields are Hq = (a), and Hq = 60.0 (b). 
Solid lines — exact solutions, open circles — approximation. 
Agreement is excellent. 



only this term determines the value of So(t). There- 
fore, we choose g\ = g% = g = Af-v3/2, and differ- 
entiate Eqs. and |S1 with respect to time once more. 
Then we see that both Eqs. El and |H1 produce the same 
result: w 7 o(t) — —woj(t) — {A\/2){wq~ i — iu 7 o), so 
that the functions w^o(t) produced by the approximate 
equations J7J and by the exact equations 10 coincide, 
provided that the initial conditions w 7 o(0) and w 7 o(0) 
also coincide. The latter condition is satisfied when all 
w a p{fy — 0, so the method described above is applicable 
only to unpolarized baths. On the other hand, for po- 
larized baths, one can use the perturbational approaches 
0,0 j so this limitation is not serious. Therefore, in or- 
der to fix the standard TDMF in case of spins 1/2, we 
just need to take Eqs.El and replace Ak by AkV3/2. It 
may seem that we just derived the standard semiclas- 
sical equations of motion, but this is not correct. For 
instance, if we take unequal spins 1 and 1/2 (So = 1 and 
Sk = 1/2) then the analysis above gives g\ — AkV3/2, 
g2 = Ak\fl /2, while the semiclassics would give g\ = 
AkV3/2, 52 = AkV2. Moreover, the initial conditions 
in our approach and in the semiclassical approach are 
different. For example, in case when the central spin 
1/2 is initially directed along the z-axis, so that the ini- 
tial density matrix p(0) = 2~ N \ |)(| | ® & =i 1, our ap- 
proach requires p({0j,(j>j},0) = (47r)~ Ar ~ 1 (l + 3cos# ), 
while the semiclassical approximation would correspond 
top({0,-,0j},O) = (4tt)- jv - 1 (5(cos6Io-1), where 6(. . .) is 
the Dirac's delta-function. 

The approach developed above gives excellent results 
at both short and long times, for different distributions of 
the coupling constants Ak, external fields Hq, and is ap- 
plicable to both longitudinal and transversal relaxation. 
A small fraction of the representative test results for a 
moderate number of bath spins (N — 15-20), is shown 
in Figs. I1I2| where we used the Chebyshev expansion 
method [14J for exact solution of the quantum problems. 
The longitudinal decay shown in Fig. ^ is typical; the 
long-time tail suggests the slow relaxation Sq(4) ~ 1/lnt, 
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Time Time 

FIG. 2: Longitudinal relaxation s§(t) of the central spin, 
couplings Ak are randomly distributed between -0.4 and 0.6, 
field Ha = 1.0. (a) central spin 1/2, N = 21 bath spins 
(b) central spin 1, iV = 19 bath spins. Solid lines — exact 
solutions, open circles — our approximation; agreement is 
excellent. Dashed lines — standard TDMF; the disagreement 
with the exact solution is significant. 
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FIG. 4: Longitudinal relaxation s§(i) of the central spin 1/2 
coupled to a bath of spins 1/2 via the anisotropic X-Y Hamil- 
tonian, external field Ho — 1.0. (a): N = 21 bath spin, Ak 
are randomly distributed between and 1.0. Solid line — 
exact solution, open circles — approximation, (b): N = 2000 
spins, the couplings Ak are the same as in Fig. but for 
smaller lattice N x = N y = 20, N z = 5. 




FIG. 3: Long-time relaxation of the central spin 1/2 cou- 
pled to a bath of 16000 spins 1/2, field Ho — 0. Graphs 
show l/so(t) as a function of Int. The coupling constants 
were calculated as At = (l/14)w(xj,), where the u(x) is the 
electron density, (a): «(x) is taken as a Gaussian with the 
half- widths d x = 8.4a, d y = 9.1a, d z = 2.2a (a is the lat- 
tice parameter), shifted from the center of the lattice by the 
vector (0.252a, 0.448a, 0.1a); (b): tt(x) is taken as an expo- 
nential function of x, with the same parameters. We used an 
extra averaging over 20 neighboring time points to decrease 
the number of realizations. The solid lines are obtained from 
raw data. 

but the results are not conclusive, since for moderate- 
sized systems s§(i — ► oo) saturates at a value far from 
zero. Thus, we have to study larger systems and to use 
the approximate method. 

We considered the systems with N — 1000-15000 spins. 
For small N = 2-5, our approximation is very crude be- 
cause the equations of motion (7J| lead to the appear- 
ance of higher-order spherical harmonics in the func- 
tion p({6j,(j)j},t) (the terms proportional to c^c^c\c v k 
etc.). These terms do not change the form of p(t) (they 
disappear after integration in Eq. 0}, but they affect 
the equations of motion for physically relevant terms, 
i.e. the actual equation of motion for uijo(t) becomes 
Wjo(t) = (A1/2)(wq 1 — Wjo) + V, where V is the con- 
tribution from the higher-order harmonics. However, the 
contribution of V into so(t) is bounded, and quickly de- 
creases for larger N, so we expect our approximation for 
srj(i) to work the better the larger N is. This is a natural 



expectation, since our method is based on TDMF, which 
works better for larger number of bath spins coupled 
to the central spin. Moreover, this assumption is well- 
confirmed by numerous numerical tests. In Fig. [21 we 
present the long-time longitudinal relaxation of an elec- 
tron spin in a model QD; we assumed that the bath spins 
1/2 are placed at the sites of a piece of a cubic lattice with 
the size N x x N y x N z , (N x = N y = 40, N z = 10, so the 
total number of bath spins N = N x N y N z = 16000). The 
long-time relaxation is extremely slow, clearly demon- 
strating the law 1/ In t (see Eil) . Our simulations show 
that the law 1/ lnf holds for unpolarized baths for differ- 
ent forms of the electron densities, i.e. for different distri- 
butions of Ak (two examples are given in Fig.|3). Our ap- 
proach also performs well (see Fig. 0J for an anisotropic 
X-Y coupling between the central and the bath spins 
H = H Q SZ + E fc A k (S$ S% +S%S y k ), which is important 
for analyzing interesting experiments of Ref. U The re- 
sults show good qualitative agreement with the exper- 
imental curves 0], but the experiments are performed 
with ~ 5-10% polarization, so that our model needs fur- 
ther development for rigorous quantitative comparison. 

Summarizing, we used the spin coherent states P- 
representation to develop a novel approach to the quan- 
tum central spin problem, based on the time-dependent 
mean field theory. The approach gives excellent agree- 
ment with the exact solutions, and is valid for a wide 
range of systems and conditions. We use it to study the 
long-time longitudinal relaxation of the electron spin in 
a quantum dot, and find that the slow decay X/hxt is 
observed in different situations. Our approach provides 
an interesting extension of the mean-field theory, and is 
applicable to many-spin central systems as well. 
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